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Abstract  In  this  investigation,  a  combined  experi¬ 
mental  and  computational  approach  with  a  Modified 
Mohr  Coulomb  (MMC)  fracture  criterion  employing 
post-initiation  element  softening  is  used  to  simulate 
stable  crack  propagation  under  Mode  I,  Mode  III 
and  combined  Mode  I/III  loading  conditions.  Results 
from  the  studies  demonstrate  that  good  correlation 
exists  between  the  measured  load-displacement  and 
the  numerically  predicted  response  when  the  stiff¬ 
ness  of  the  specimen  fixture  is  included  in  the  FE 
model.  The  numerical  results  were  able  to  capture 
most  of  the  experimentally  observed  features  dur¬ 
ing  crack  propagation,  such  as  through-thickness  slant 
fracture,  necking,  tunneling  and  local  specimen  twist, 
thus  confirming  that  the  MMC  criterion  is  suitable  for 
predicting  in-plane  and  out-of-plane  tearing  of  sheets. 
It  was  found  that  in  order  to  predict  correctly  the 
load-displacement  curve  as  well  as  the  fracture  plane, 
different  amount  of  softening  is  needed  for  Mode  I  and 
Mode  III  loading  cases.  This  observation  can  be  jus¬ 
tified  on  the  micro-mechanical  level,  while  there  is  a 
competition  between  the  mechanisms  of  dimple  and 
shear  fracture. 
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1  Introduction 

The  problem  of  ductile  crack  initiation  and  propagation 
for  a  pre-cracked  specimen,  such  as  Compact  Tension 
(CT)  specimens,  subjected  to  a  single  mode  of  load¬ 
ing  (e.g.,  Mode  I)  has  received  a  great  deal  of  attention 
in  the  literature.  Several  methods,  such  as,  cohesive 
elements,  CTOD,  CTOA,  and  the  J-integral  can  pre¬ 
dict  fracture  response  with  varying  degrees  of  accu¬ 
racy.  Even  though  this  area  of  fracture  mechanics  is 
well  established,  there  is  much  less  information  in  the 
literature  regarding  combined  Mode  I/II  or  Mode  I/III 
response. 

Parallel  advances  were  made  in  the  experimental 
part  of  this  research.  Two  types  of  clamping  fixtures 
developed  for  Mode  I/II  loading  and  Mode  I/III  load¬ 
ing,  suitable  for  standard  tensile  rigs  are  shown  in  Fig.  1 . 

Presently,  the  displacement  and  the  strain  fields  are 
measured  by  digital  correlation  system.  Tearing  tests  on 
pre-crack  specimens  of  the  type  shown  in  Fig.  1  provide 
a  wealth  of  information  needed  to  develop  a  suitable 
fracture  model.  Another  type  of  equipment  suitable  for 
studying  the  mixed  Mode  I/II  fracture  is  the  modified 
Arcan  fixture  (Amstutz  et  al.  1995, 1997).  In  all  of  these 
studies,  tests  were  performed  on  compact  tension  spec¬ 
imens  so  that  crack  initiation  was  not  part  of  either  the 
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Fig.  1  Special  fixtures  for 
combine  Mode  I/II  (left, 
Aoki  et  al.  1990)  and  Mode 
I/III  (right,  Yan  et  al.  2009a) 
loading 


experiments  or  the  predictions.  Based  on  the  experi¬ 
mental  measurements,  which  included  both  the  mixed 
mode  crack  opening  displacement  vector,  COD,  and 
the  surface  strain  fields,  a  mode  transition  from  tension- 
dominated  to  shear-dominated  was  identified.  Based  on 
these  observations,  a  mixed  mode  COD  criterion  was 
proposed  to  predict  both  the  initial  crack  growth  and 
direction  of  crack  extension  (Ma  et  al.  1999;  Sutton 
et  al.  2000a, b).  As  part  of  these  studies,  the  investi¬ 
gators  noted  that  components  of  COD  are  most  likely 
a  function  of  constraint  at  the  crack  tip,  resulting  in  a 
series  of  articles  that  focused  on  the  underlying  physics 
related  to  the  combined  effect  of  constraint  and  effec¬ 
tive  stress  (Sutton  et  al.  2002;  Zuo  et  al.  2004,  2005). 
This  work  was  later  extended  into  experimental,  theo¬ 
retical  and  numerical  studies  to  describe  the  effect  of 
tension-torsion  loading  on  fracture  of  thin  sheet  mate¬ 
rials  (Sutton  et  al.  2001,  2007;  Lan  et  al.  2007;  Yan 
et  al.  2007, 2009a, b)  with  additional  focus  on  slant  frac¬ 
ture  and  whether  its  presence/absence  can  be  predicted 
using  macroscopic  quantities  such  as  effective  plastic 
strain  or  effective  stress  (Mahgoub  et  al.  2003 ;  Lan  et  al. 
2006). 

In  several  recent  publications,  comparison  was  made 
between  the  numerical  prediction  and  experimental 
results.  Sutton  et  al.  (2000a, b)  employed  their  gener¬ 
alized  CTOD  method  to  predict  crack  growth  during 
mixed  Mode  I/II  fracture,  with  the  direction  of  crack 
growth  predicted  using  a  specified  ratio  of  COD  com¬ 
ponents  and  the  onset  of  crack  growth  defined  by  the 
magnitude  of  the  critical  COD. 


The  Cohesive  Zone  Model  (CZM)  provides  also  a 
practical  way  to  predict  mixed-mode  fracture  (Ortiz  and 
Pandolfi  1999).  De-Andre  et  al.  (1999)  demonstrated 
an  application  of  cohesive  elements  to  fatigue  crack 
growth  using  a  three-dimensional  cohesive  element 
with  irreversible  cohesive  laws  to  track  three- 
dimensional  fatigue  crack  fronts.  The  simulation  results 
were  compared  with  the  axial  fatigue  tests  on  alu¬ 
minum  shafts  performed  by  Thompson  and  Sheppard 
(1992).  Simonsen  and  Tornqvist  (2004),  reported  on 
a  combined  experimental-numerical  procedure  to  cali¬ 
brate  macroscopic  crack  propagation  criteria  (Rice  and 
Tracey  (1969)  and  a  constant  equivalent  strain  crite¬ 
rion)  for  large-scale  structures,  with  an  application  to 
ship  grounding.  Oliver  (1996a, b),  Teng  (2008)  and  Xue 
(2007a)  reported  on  successful  application  of  contin¬ 
uum  damage  mechanics  in  simulating  initiation  and 
propagation  of  cracks.  An  alternative  approach  is  being 
developed  by  Belytschko  and  his  followers  using  the 
concept  of  Extended  Finite  Element  method  (XFEM) 
(Moes  and  Belytschko  2002;  Gravouil  et  al.  2002). 
Other  less  general  numerical  methods  are  also  pro¬ 
posed,  for  example,  (Morrissey  and  Geubelle  1997; 
Gao  and  Klein  1998;  Oliver  1996a, b). 

More  recently,  Bai  and  Wierzbicki  (2010)  extended 
the  Mohr-Coulomb  failure  criterion  (MMC)  and  intro¬ 
duced  an  explicit  dependence  of  the  equivalent  strain 
to  fracture  on  both  the  stress  triaxiality  and  Lode  angle 
parameter.  The  Lode  angle  6  is  proportional  to  the  nor¬ 
malized  third  invariant  of  stress  deviator.  This  fracture 
criterion  includes  three  parameters  which  are  calibrated 
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from  a  limited  number  (minimum  three)  of  experi¬ 
ments.  The  MMC  method  is  an  alternative  to  the  widely 
used  GTN  model,  which  requires  different  types  of 
parameters  to  be  determined  from  experiments. 

Several  people  in  the  ICL  lab  contributed  to  the 
development  of  the  MMC  fracture  model,  which  is  the 
current  working  model  in  the  lab.  Bao  and  Wierzbicki 
(2004)  constructed  an  empirical  fracture  diagram  in  a 
wide  range  of  stress  triaxiality  based  on  the  results 
of  36  meticulous  experiments  and  numerical  simula¬ 
tions  of  several  types  of  specimens.  More  recently,  Xue 
(2007b)  and  Xue  and  Wierzbicki  (2008)  introduced  the 
dependence  of  the  fracture  strain  on  the  Lode  angle  and 
constructed  a  general  form  of  an  asymmetric  fracture 
model.  Bai  and  Wierzbicki  (2010)  showed  that  the  three 
dimensional  fracture  locus  could  be  derived  in  a  sim¬ 
ple  and  elegant  way  from  the  modified  Mohr-Coulomb 
failure  criterion.  The  post-failure  model  was  developed 
and  applied  to  the  plane  strain  initiation  and  propa¬ 
gation  of  cracks  by  and  Li  and  Wierzbicki  (2010a). 
Finally,  MMC  model  was  calibrated  for  several  grades 
of  Advanced  High  Strength  Steels  (AHSS).  The  valida¬ 
tion  of  this  model  in  sheet  metal  forming  operations  was 
presented  in  three  recent  publications  (Li  et  al.  2010b; 
Luo  and  Wierzbicki  2010).  For  the  sake  of  economy 
space,  the  details  of  the  derivation  of  MMC  are  not 
given  in  this  paper.  The  interested  reader  is  referred  to 
the  above  original  publications. 

It  is  noted  that  all  plasticity  and  fracture  parameters 
used  in  the  numerical  simulations  were  obtained  from 
an  entirely  independent  set  of  experiments  performed 
at  MIT  (Beese  et  al.  2010),  whereas  the  experiments 
being  simulated  were  performed  at  USC  and  NASA. 
Thus,  the  comparison  of  the  predicted  and  measured 
load-displacement  curve  is  an  essential  part  of  the  val¬ 
idation  process  for  the  simulations. 

In  this  study,  the  MMC  stress-strain  based  fracture 
criteria  is  used  to  simulate  crack  propagation  under 
nominally  Mode  I,  Mode  III  and  combined  Mode  I/III 
conditions.  An  extensive  comparison  is  made  with 
experiments  performed  by  Yan  et  al.  (2009a, b)  using 
A16061-T6  sheet  specimens. 


2  Material  model 


Fig.  2  Stress-strain  data  for  A16061-T6  (all  simulations  in  this 
paper  used  the  averaged  MIT  data) 


ogy  (MIT)  and  the  University  of  South  Carolina  (USC). 
Though  the  material  used  at  both  laboratories  came 
from  different  batches,  the  measured  hardening  curves 
(red  at  MIT,  blue  at  USC)  are  quite  similar,  as  shown 
in  Fig.  2.  Baseline  experiments  performed  at  MIT  by 
extracting  dog-bone  specimens  at  0°,  45°  and  90°  from 
the  sheet  rolling  direction  indicated  the  presence  of 
weak  anisotropy.  In  this  study,  all  simulations  were 
performed  using  a  standard  J2  plasticity  model  with 
the  averaged  MIT  stress-strain  curve  and  without  con¬ 
sidering  anisotropy  in  the  stress-strain  response.  The 
finite  element  calculations  were  performed  by  ABA- 
QUS 6. 8/Explicit  which  is  a  nonlinear  FE  code  for  large 
plastic  strain. 


3  Fracture  model 

The  fracture  model  used  in  these  simulations  was 
described  in  previous  work  (Bai  and  Wierzbicki  2010; 
Beese  et  al.  2010;  Li  and  Wierzbicki  2010a;  Li  et  al. 
2010b;  Li  and  Wierzbicki  2009).  According  to  this 
model,  fracture  in  a  given  element  is  assumed  to  occur 
when  the  accumulated  normalized  plastic  strain  reaches 
unity 


D  = 


d£r 


(la) 


All  specimens  were  manufactured  using  6061-T6  alu¬ 
minum  alloy,  which  was  extensively  studied  by  the 
present  authors  at  Massachusetts  Institute  of  Technol- 


where  sp  is  the  equivalent  plastic  strain,  r]  =  ^  stands 
for  stress  triaxiality,  am  is  the  mean  stress,  a  is  the 
equivalent  stress,  and  0  is  the  Lode  angle  parameter, 
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which  is  defined  as 

0  =  1  —  —  arccos  ((=)3)  (lb) 

where  r  is  the  third  invariant  of  the  deviatoric  stress 
tensor.  The  meaning  of  the  denominator  of  the  inte¬ 
grand  of  Eq.  1  becomes  clear  when  we  consider  a  class 
of  proportional  loading  paths  for  which  the  parame¬ 
ters  77  and  0are  constant.  In  this  case,  Eq.  1  is  reduced 
to  Sf  =  f(rj,6)  after  integration.  Thus  the  normal¬ 
izing  function  in  the  denominator  of  Eq.  1  represents 
the  magnitude  of  the  equivalent  fracture  strain  which 
is  reached  by  all  possible  proportional  loading  paths. 

The  general  expression  for  the  equivalent  strain  to 
fracture  as  a  function  of  the  stress  triaxiality  and  Lode 
angle  parameters,  derived  from  the  MMC  model  is 


Sf(rj,6) 


There  are  three  fracture  constants  in  the  above  equation, 
Ci  C2  and  C3that  could  be  determined  from  a  mini¬ 
mum  of  three  tests.  In  the  present  problem,  five  differ¬ 
ent  tests  were  performed  which  increases  the  accuracy 
of  the  fracture  locus. 

Since  the  specimen  used  in  the  Mode  I/III  exper¬ 
iments  is  relatively  thin,  plane  stress  is  a  reasonable 
assumption.  In  this  case,  the  MMC  fracture  model 
yields  the  following  set  of  equations  (see  Bai  and 
Wierzbicki  2010). 


Sf(v)  = 


■fi 


+C 1 
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(3) 


/l  =  COS 
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27  / 

2  All 
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2  i\ii 
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Table  1  Calibrated  material  and  fracture  parameters  for 
A16061-T6 


A 

n 

Ci 

C2 

C3 

438  MPa 

0.07 

0.06 

228  MPa 

0.93 

Fracture  Initiation  Crack  propagation  Element  failure 

Fig.  3  The  concept  of  the  simulation  of  post-initiation  behavior 
of  an  element  (Li  and  Wierzbicki  2010a) 


In  the  above  equations,  n  is  the  hardening  exponent 
obtained  using  the  uniaxial  data  in  Fig.  2,  The  ampli¬ 
tude  of  the  hardening  law  is  denoted  by  A  and  the  three 
model  parameters  Ci,  C2,  C3  related  to  the  functional 
form  of  the  critical  fracture  strain  are  to  be  determined 
from  a  series  of  fracture  experiments.  A  set  of  plastic¬ 
ity  and  fracture  parameters,  determined  by  Beese  et  al. 
(2010)  for  A16061-T6  and  used  in  the  present  simula¬ 
tions,  are  shown  in  Table  1. 

The  graphic  representation  of  Eq.  3  is  shown  in 
Fig.  19.  There  are  four  segments  of  the  fracture  locus, 
each  represented  by  a  different  color  in  Fig.  19. 


4  Post-initiation  behavior  of  an  element 

In  order  to  describe  the  formation  of  slant  fracture, 
some  type  of  post-initiation  softening  must  be  assumed 
(Li  and  Wierzbicki  2010a).  One  can  define  two  types 
of  softening,  one  originating  from  material  weakening 
due  to  void  growth  and  linkage.  This  type  of  softening 
has  been  extensively  studied  in  the  literature.  In  addi¬ 
tion  to  the  existence  of  softening  process  at  the  material 
level,  one  can  interpret  the  softening  at  the  level  of  a 
finite  element.  This  is  illustrated  in  Fig.  3. 

Fracture  is  assumed  to  initiate  at  one  point  within 
an  element.  Because  of  a  finite  size  of  the  element, 
additional  work  is  needed  to  propagate  the  crack 
through  the  element.  The  local  material  (finite  ele¬ 
ment)  gradually  loses  strength  as  crack  propagates 
through  the  element  and  reaches  the  element  bound¬ 
ary  when  total  separation  finally  occurs.  The  above 
interpretation  of  “structural”  softening  has  led  to  a 
realistic  model  of  the  fracture  process  that  allowed 
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the  investigators  to  successfully  predict  the  onset  of 
slant  fracture  (Li  and  Wierzbicki  2010a).  In  the  present 
paper,  the  post-initiation  model  available  in  ABAQUS 
V6.8  will  be  used.  This  model  involves  only  one  scalar 
parameter  u  /  that  may  change  from  one  loading  case  to 
the  other.  For  example,  we  have  found  that  larger  values 
of  u  f  should  be  used  for  Mode  I  and  smaller  values  for 
Mode  III  fracture,  with  instantaneous  fracture  corre¬ 
sponding  to  uf  =  0.  The  larger  the  parameter  uf,  the 
longer  post-initiation  lasts.  In  our  present  numerical 
simulations,  different  values  of  Uf  were  used,  ranging 
from  le~5  to  0.1  (see  Sect.  10.2). 

5  Description  of  experiments 

All  mixed  mode  I/III  experiments  were  performed  in  an 
MTS  810  load  frame  under  displacement  control,  see 
Fig.  4.  The  experimental  set-up  consists  of  securely 
attaching  the  specimen  to  two  quarter  circle  fixtures. 
The  fixtures  are  connected  to  the  loading  frame  through 
pins  that  are  inserted  in  loading  clevises,  see  Fig.  5.  The 
clevises  in  turn  are  clamped  into  the  hydraulic  grips  of 
the  loading  frame.  Depending  on  the  position  of  the 
pin  with  respect  to  the  fixtures,  one  can  generate  Mode 
I,  Mode  III  and  combined  Mode  I/III  far  field  loading 
conditions.  The  position  of  the  pins  is  characterized 
by  the  loading  angle  O,  defined  in  Fig.  4.  The  value 
O  =  0o  corresponds  to  Mode  I  loading,  while  O  =  90° 
provides  nominally  far-field  Mode  III  conditions.  The 
intermediate  angles  O  =  30°  and  =  60°  generate 
combined  Mode  I/III  fracture. 

In  this  work,  three-dimensional  Digital  Image  Cor¬ 
relation  (3D-DIC)  is  used  to  obtain  surface  displace¬ 
ment  measurements  in  the  crack  tip  region  in  these 
experiments.  The  3D-DIC  system  is  shown  in  Fig.  4. 
Before  performing  each  experiment,  a  high  contrast 
random  pattern  of  speckles  is  applied  to  each  speci¬ 
men  and  the  3D-DIC  system  is  positioned  to  view  the 
crack  tip  region.  Prior  to  performing  each  experiment, 
both  cameras  are  synchronized  and  calibrated.  A  photo¬ 
graph  of  the  as-painted  undeformed  specimen,  as  well 
as  a  photograph  of  a  deformed  and  fracturing  specimen 
with  O  =  60°,  are  shown  in  Fig.  6. 

6  FE  modeling 

Two  models  were  developed  by  discretizing  the 
deformable  part  of  the  specimen,  using  either  solid 


1 — Camera  0,  2 — Camera  1T  3 — Translation  Stage,  4 — Hydraulic 
Grip  5 — Clevis,  6 — Fixture,  7 — Specimen,  8 — Light  Source,  9- 
Rotation  Constraint,  10 — extension  bar 


Fig.  4  Experimental  setup  for  combined  tension-torsion  loading 
of  thin  specimen 


Fig.  5  Photos  of  specimen  fixture  (left)  and  the  clevis  connecting 
the  fixture  with  the  cross-head  of  the  testing  machine  (right) 

or  shell  elements.  Both  models  have  a  refined  cen¬ 
tral  part  along  the  expected  path  of  the  crack.  Out¬ 
side  this  region,  a  much  coarser  mesh  was  used.  The 
solid  element  is  eight  nodes  with  reduced  integration 
(C3D8R  in  ABAQUS).  The  shell  element  is  four  nodes 
with  reduce  integration  (S4R  in  ABAQUS).  In  the 
base  FE  model,  the  solid  element  size  was  very  small, 
equal  to  0.0625  mm  in  the  critical  region  along  the 
crack  path.  Other  models  with  larger  mesh  size  will 
be  discussed  in  the  following  sections.  In  all  cases 
there  were  5  elements  through  the  thickness.  Alto¬ 
gether  there  were  about  38,000  elements  in  the  base 
model.  A  representative  undeformed  mesh  is  shown  in 
Fig.  7.  It  is  remarked  that  the  choice  of  the  number  of 
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Fig.  6  Dimensions  of  the 
undeformed  specimen  (left) 
and  the  deformed  specimen 
with  about  20  mm  crack 
extension  during  a  mixed 
mode  I/III  experiment  at 
loading  angle  <J>  =  60° 


elements  through  the  thickness  was  made  as  a  com¬ 
promise  between  accuracy  and  CPU  time.  Five  ele¬ 
ments  through  the  thickness  are  sufficient  to  simulate 
slant  fracture,  which  was  one  of  the  purposes  of  this 
paper. 

For  modeling  purposes,  the  parts  of  the  specimen 
bolted  to  the  lower  and  upper  fixtures  were  assumed 
to  be  “rigid”  and  represented  by  two  reference  points 
(see  blue  dots  in  Fig.  7).  The  two  reference  points  are 
tied  to  the  interface  of  the  rigid  part  and  the  deform¬ 
able  part  of  the  specimen.  Since  the  lower  fixture  is 
much  stiffer  than  the  specimen,  the  investigators  opted 
to  use  a  ‘U -joint  connector’  option  available  in  ABA- 
QUS  to  conveniently  represent  the  behavior  of  the  con¬ 
nection  to  the  clevises,  see  Fig.  8.  It  should  be  noted 
that  a  major  advantage  of  using  the  connector  concept 
to  describe  the  overall  behavior  of  the  lower  fixture  is 
that  ABAQUS  allows  the  investigator  to  incorporate 
“stiffness  variables”  (e.g.  rotational  stiffnesses  Ki  and 
K2).  These  variables  provide  the  opportunity  to  sim¬ 
plify  the  incorporation  of  complex  joint/fixture  behav¬ 
ior  into  the  simulations.  For  example,  the  investigators 
observed  effects  such  as  slippage  and  small  specimen 
rotation  within  the  “rigid”  bolted  connection  that  would 
affect  the  overall  response  but  are  difficult  to  model 
accurately  unless  model  parameters  are  used  to  repre¬ 
sent  these  effects. 

The  control  variables  in  the  experiment  and  simu¬ 
lation  are  the  vertical  displacements  of  the  pins.  The 
output  parameters  of  the  finite  element  simulation  are 
the  relationship  between  the  total  reaction  force  and 
the  pin  displacement,  the  trajectory  of  the  crack  and  a 
complete  set  of  information  about  the  stress,  strain  and 
displacement  field  around  the  crack  tip.  This  allows  for 


a  thorough  comparison  with  experiments  where  many 
of  the  above  parameters  were  measured. 

7  Formulation  of  boundary  conditions 

As  shown  in  Fig.  4,  the  thin  specimen  is  firmly  clamped 
to  the  fixture  (component  6  in  Fig.  4)  with  eight  bolts, 
the  clevis  (component  5)  is  attached  to  the  hydraulic 
grip  (component  4)  via  high  pressure  clamping  plates 
and  the  upper  fixtures  (component  6)  is  attached  to 
the  upper  clevis  through  an  “extension  bar”  (compo¬ 
nent  10).  Both  experimental  observations  and  inspec¬ 
tion  of  the  design  drawings  suggest  that  there  is  a 
fair  amount  of  flexibility  in  the  overall  experimental 
configuration  (e.g.,  pin  connections  between  the  clevis- 
extension  bar-upper  lower  fixture;  slight  specimen  rota¬ 
tion  due  to  slippage  within  lower  fixture).  This  sug¬ 
gests  that,  while  the  displacements  are  continuous  at 
all  junctures,  finite  rotational  flexibility  must  be  con¬ 
sidered  to  account  for  this  effect  when  modeling  the 
overall  configuration.  In  fact,  our  simulations  clearly 
show  that  the  load-displacement  curves  and  displace¬ 
ment-crack  propagation  results  are  highly  sensitive  to 
the  rotational  boundary  conditions,  requiring  inclusion 
of  the  deformable  connections  into  the  computational 
model. 

Figure  8  shows  the  FE  model  for  the  combined  load¬ 
ing  cases  with  loading  angle  O.  Rotational  flexibility 
of  the  clevis  is  modeled  via  the  U-joint  connectors 
between  the  pins  and  the  reference  points  of  the  rigid 
part  of  the  specimen,  using  the  stiffness  values  defined 
in  the  Appendix  and  constraining  the  translational 
degrees  of  freedom. 
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Fig.  7  Mixed  mode  I/III 
fracture  specimen  and  our 
finite  element  model 


Fig.  8  Specimen  and 
loading  configuration  for 
combined  Mode  I/III 
loading 


Using  the  global  ( laboratory )  coordinate  system 
shown  on  the  clevis  in  Fig.  9,  the  bending  moment  vec¬ 
tor,  M,  has  components  in  the  global  ( laboratory )  sys¬ 
tem  (Mi ,  M2,  M3).  If  all  rotations  are  constrained,  then 
all  three  components  of  the  bending  moment  would  be 
transmitted  to  the  specimen.  In  addition  to  the  global 
coordinate  system,  a  local  ( material )  coordinate  sys¬ 


tem  is  defined  which  is  attached  to  the  specimen,  as 
shown  in  Fig.  9.  In  this  coordinate  system,  the  bending 
moment  vector  is  M[_  with  components  (M[ ,  M^,  M3). 
Assuming  zero  friction  between  the  pin  and  the  cle¬ 
vis,  then  M3  =  M3  ~  0  and  the  3-D  problem  is 
reduced  to  the  2-D  planar  transformation  of  bending 
moments. 
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Fig.  9  Definition  of  global  and  local  coordinate  systems  a 
Global  (laboratory)  coordinate  system,  b  local  (material)  coor¬ 
dinate  system 
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Letting  K'2  =  0  (free  clevis  pin  rotation),  the  above 
system  reduces  to 

m[  =  k[o[ 

M2  =  K'20f2  (7) 

In  order  to  study  how  sensitive  the  load-displacement 
curves  are  to  the  stiffness  matrix,  limiting  cases  of  the 
springs,  pinned  or  fixed  (zero  stiffness  or  infinite  stiff¬ 
ness)  are  investigated  for  both  nominally  Mode  I  and 
Mode  III  loading.  The  actual  values  of  K[  and  K2  are 
between  the  limiting  cases  and  are  estimated  in  the 
Appendix. 


M[  =  M\  cos  +  M2  sin  O 

M2  =  —M\  sin  O  +  M2  cos  O  (4) 

where  the  relative  position  of  the  two  coordinate  sys¬ 
tems  is  uniquely  defined  by  the  angle  O  shown  in  Figs.  4 
and  8. 

The  moment  M[  represents  an  in-plane  bending 
moment  in  the  specimen  whereas  M2  represents  the 
expected  out-of-plane  bending.  It  is  noted  that  one 
physical  source  for  the  moment  M2  is  the  resistance 
of  the  arms  of  the  fork  to  twist  (Appendix).  Of  partic¬ 
ular  concern  is  the  component  Mi  which,  when  trans¬ 
ferred  to  the  specimen  system,  contributes  to  ‘in-plane 
bending’  of  the  specimen.  The  component  Mi  could 
result  from  any  combination  of  (a)  incremental  slippage 
within  the  bolted  fixture- specimen  connection  causing 
in-plane  rotation  of  the  separate  arms  of  the  specimen, 
(b)  non-zero  clearance  between  lower  fixture  and  cle¬ 
vis  allowing  rigid  body  rotation  of  fixture  (see  Fig.  9) 
and/or  (c)  finite  flexural  stiffness  of  each  clevis  fork 
arm. 

It  is  assumed  that  the  moment  vector  M_'  is  related 
to  the  corresponding  rotation  0'  through  a  linear  trans¬ 
formation: 

M  =  K0_  (5a) 

Mr  =  K[  0',  (5b) 

where  K_  and  K[  are  stiffness  matrices  respectively  in 
the  global  and  material  coordinate  systems.  Assum¬ 
ing  uncoupled  rotational  stiffness  factors,  the  stiffness 
matrices  K_  and  K*_  become  diagonal: 


8  Finite  element  simulations 

Since  MMC  is  a  strain  based  fracture  criteria,  the  col¬ 
ors  in  all  of  the  contours  of  the  FE  results  shown  in 
this  paper  represent  the  equivalent  plastic  strain,  red  is 
maximum  and  blue  is  minimum. 

8.1  Effects  of  stiffness  parameters 

Four  simulations  of  Mode  I  tensile  tests  were  per¬ 
formed  using  different  values  for  the  constant  Ki  using 
the  eight-node  solid  element  with  reduce  integration 
(C3D8R  in  ABAQUS)  and  the  standard  mesh  with  ~ 
38,000  elements.  The  corresponding  load-displacement 
curves  and  the  experimental  data  are  shown  in  Fig.  10a. 
The  line  corresponding  to  the  infinite  stiffness  (full 
rotational  fixity)  over-predicts  the  test  results  by  a  fac¬ 
tor  of  two.  At  the  same  time,  by  releasing  the  rotational 
degree  of  freedom  (K\  =  0)  the  specimen  resistance  is 
underestimated.  The  simulations  of  the  two  extreme 
values  of  the  planar  rotational  stiffness  demonstrate 
that  the  clamping  device  of  the  specimen  can  not  be 
assumed  to  be  perfectly  rigid.  Additional  simulations 
for  two  intermediate  values  of  the  rotational  constant 
K 1  introduced  a  dramatic  improvement  in  the  accuracy 
of  the  numerical  prediction. 

It  is  noted  that  the  Mode  I  response  is  practi¬ 
cally  a  planar  problem  so  that  torsional  stiffness  of 
the  clevis  K2  does  not  influence  the  solution.  By 
contrast,  the  torsional  stiffness  has  an  effect  on  all 
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Fig.  10  The  effects  of  finite 
stiffness  a  Ki  effect  for 
tension  loading;  b  K2  effect 
for  torsional  loading;  c  K2 
effect  for  O  =  30°  loading; 
d  K2  effect  for  <t>  =  60° 
loading 
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other  loading  cases  (see  Fig.  lOb-d).  For  example,  the 
response  for  shear  loading  was  simulated  using  the 
estimated  value  of  the  bending  stiffness  for  the  fix¬ 
ture  K\  =  2,000  Nm  (see  the  Appendix)  and  three 
different  values  of  the  torsional  spring  constants  (see 
legend  of  Fig.  10b).  Generally,  the  curve  correspond¬ 
ing  to  the  estimated  value  of  the  torsional  stiffness 
K2  =  800  Nm,  lies  between  the  two  extreme  cases 
of  zero  and  infinite  stiffness.  Clearly  the  effect  of  tor¬ 
sional  stiffness  is  pronounced  but  is  not  as  strong  as 
the  effect  of  K\  on  Mode  I  response.  Similar  conclu¬ 
sions  can  be  drawn  for  the  combined  loading  cases 
shown  in  Fig.  10c  and  d.  Based  on  the  above  para¬ 
metric  study,  the  values  K\  =  2,000,  K2  =  800  Nm 
and  K3  =  0  will  be  used  for  other  loading  cases  (30 
and  60  degree). 


8.2  Comparison  with  experimental  measurements 

All  simulations  were  performed  with  ABAQUS6.8 
using  the  eight-node  solid  element  C3D8R  and  a  stan¬ 
dard  mesh  (see  Fig.  7). 


8.2.1  Load- displacement 

Comparison  between  the  measured  and  calculated 
load-displacement  curve  for  all  four  cases  is  shown 
in  Fig.  11,  taking  the  optimum  values  of  the  stiffness 
matrix  obtained  from  the  above  parametric  study.  The 
correlation  is  good  for  all  loading  angles,  providing 
confidence  that  the  simulations  are  capable  of  ade¬ 
quately  transferring  global  load  conditions  to  the  spec¬ 
imen  for  a  wide  range  of  loading  angles,  initiation  of 
crack  growth  and  a  large  range  of  crack  extensions. 
The  inserted  color  coded  pictures  in  Fig.  1 1  show  the 
deformed  shape  of  the  specimen  in  an  intermediate 
stage  of  crack  propagation  for  all  four  loading  angles. 

8.2.2  Crack  path 

Within  the  scope  of  the  present  element  deletion  tech¬ 
nique,  crack  growth  is  considered  as  a  sequence  of 
crack  initiation.  In  order  to  delete  a  given  element,  the 
accumulated  damage  should  reach  the  critical  value 
and  the  load  carrying  capacity  must  vanish  in  the 
post-initiate  range.  It  should  be  noted  that  in  the  above 
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Fig.  11  Results  of 
numerical  simulations  for 
four  values  of  the  loading 
angle  O 
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approach,  the  constraints  at  the  crack  tip  (defined  by 
stress  triaxiality  and  Lode  angle  parameter  may  change 
in  the  course  of  the  deformation  of  a  given  element.  The 
predicted  crack  path  is  straight  because  of  the  assumed 
symmetry  of  the  boundary  conditions  on  both  sides  of 
the  specimens  and  in  generally  good  agreement  with 
experimental  evidence.  It  must  be  noted  that  for  mode 
III  loading  case,  Yan  et  al.  (2009a)  observed  a  slightly 
curvilinear  crack  growth  for  O  =  90° .  In  their  exper¬ 
iment,  an  extension  bar  was  added  to  the  load  train 
between  the  upper  clevis  and  fixture  to  provide  better 
visual  access  for  the  stereo  vision  measurements.  Such 
a  modification  of  the  gripping  system  introduced  asym¬ 
metry  to  the  specimen  loading  by  lowering  the  bend¬ 
ing  stiffness  and  changing  the  torsional  stiffness  of  the 
upper  clevis.  Since  it  is  difficult  to  quantify  the  amount 
of  asymmetry,  the  authors  opted  for  maintaining  sym¬ 
metry  in  the  simulations  by  assuming  that  the  lower  and 
upper  portions  of  the  load  train  have  the  same  stiffness 
values  for  all  loading  angles. 

9  Comparison  between  solid  and  shell  element 
model  predictions 

For  large  structures,  the  solid  element  models  used  in 
the  previous  sections  are  computationally  expensive. 


To  determine  whether  the  tearing  process  can  be  mod¬ 
eled  with  reasonable  accuracy  by  shell  elements,  the 
present  specimen  was  modeled  by  four-node  shell  ele¬ 
ments  with  reduced  integration  (S4R).  Five  integration 
points  were  taken  through  the  wall  thickness.  All  other 
input  parameters  (e.g.,  planar  mesh  size,  post-initiation 
parameter  u  /,  and  material  data)  were  taken  to  be  the 
same  as  used  with  the  solid  element  model.  The  ele¬ 
ment  deletion  rule  for  shell  element  was  implemented, 
such  that  an  element  is  deleted  when  all  the  integration 
points  lose  their  strength. 

Figures  12  and  1 3  presents  a  comparison  of  the  mod¬ 
eling  results  for  Mode  I  and  Mode  III  loading,  respec¬ 
tively.  As  shown  in  Figs.  12  and  13,  both  solid  and 
shell  models  adequately  capture  the  load-displacement 
behavior  before  fracture  initiation.  However,  the  shell 
element  model  do  not  accurately  model  crack  propaga¬ 
tion  for  either  Mode  I  or  Mode  III  loading.  For  Mode 
I  loading,  the  shell  element  model  predicts  that  crack 
growth  starts  much  earlier  than  the  solid  element  model; 
the  load  reached  only  half  of  the  experimentally  mea¬ 
sured  peak  load,  and  the  process  of  Mode  I  crack  prop¬ 
agation  is  unstable  after  growth  initiates.  For  the  Mode 
III  loading  case,  the  shell  element  model  correctly  pre¬ 
dicted  crack  initiation  and  the  corresponding  peak  force. 
However,  large  discrepancies  appear  later  during  the 
crack  propagation  phase,  as  evidenced  in  Fig.  13. 
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Fig.  12  Simulation  results  of  shell  element  model  and  solid  ele¬ 
ment  model  for  Mode  I  loading 
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Fig.  13  Simulation  results  of  shell  element  model  and  solid  ele¬ 
ment  model  for  Mode  III  loading 

Based  on  our  preliminary  studies,  it  is  clear  that  the 
shell  element  models  are  not  suitable  for  predicting 
tearing  process  in  sheets  and  hence  all  further  model 
predictions  will  be  based  on  the  solid  element  model. 

10  Parametric  study 

In  this  section,  sensitivity  analysis  is  performed  on  the 
effect  of  several  modeling  parameters  on  the  finite  ele¬ 
ment  solution. 

10.1  Influence  of  mesh  size 

Several  solid  model  simulations  were  performed  while 
maintaining  a  constant  size  of  the  solid  element  through 
the  thickness  and  along  the  crack  path  equals  0.4  mm. 
At  the  same  time,  three  mesh  sizes  were  introduced 


in  the  vertical  direction  equal  respectively  to  0.0625, 
0.25  and  0.5  mm.  For  the  case  of  Mode  I  fracture,  there 
was  little  effect  of  mesh  size  on  the  load-displacement 
results.  For  Mode  III  response,  there  is  a  somewhat 
larger  difference  in  specimen  response  between  those 
three  sizes  of  meshes,  see  Fig.  14. 

In  addition,  there  is  a  clear  convergence  of  the  results 
of  the  FE  simulation  for  Mode  I  loading  with  diminish¬ 
ing  mesh  size.  In  the  case  of  Mode  III  loading,  the  con¬ 
vergence  is  much  slower.  It  can  be  concluded  that  the 
effect  of  mesh  size  of  the  global  structure  response  with 
the  present  model  is  small.  At  the  same  time,  the  mesh 
size  does  influence  the  local  response  and  in  particu¬ 
lar,  the  crack  orientation.  The  experimentally  observed 
slant  fracture  behavior  cannot  be  predicted  when  the 
mesh  size  is  too  coarse.  This  is  illustrated  in  Fig.  15, 
which  compares  results  from  two  mesh  sizes  in  the 
vertical  direction,  0.25  and  0.0625  mm.  In  the  case  of 
0.25  mm  mesh,  the  fracture  surface  was  flat  and  normal 
to  the  plane  of  the  specimen,  as  shown  in  Fig.  15a.  In 
order  to  generate  slant  fracture,  a  much  smaller  verti¬ 
cal  mesh  size  (0.0625  mm)  is  necessary,  as  shown  in 
Fig.  15b. 

10.2  Influence  of  post-initiation  parameter  u / 

After  crack  growth  occurs,  elements  near  the  initial 
growth  site  begin  to  gradually  lose  strength,  resulting 
in  “structural”  softening.  It  was  shown  by  Xue  (2007a) 
and  Li  and  Wierzbicki  (2010a),  that  the  presence  of 
softening  is  necessary  to  simulate  slant  fracture.  Two 
extreme  values  of  the  post-initiation  parameter  u  /  were 
assumed  to  show  the  sensitivity  of  the  solution;  Uf  = 
0.1  mm  and  uf  =  1^  —  5  mm.  As  shown  in  Fig.  16, 
changing  Uf  by  four  orders  of  magnitude  produces  only 
a  moderate  effect  on  the  load  displacement  curve.  It  was 
found  that  a  larger  value  is  more  appropriate  for  Mode 
I  response,  while  a  smaller  value  is  needed  for  Mode 
III  simulations.  Indeed,  the  best  fit  of  the  experimental 
curve  is  obtained  using  uf  =0.1  mm  for  Mode  I,  and 
Uf  =  0.01  mm  for  Mode  III. 

The  magnitude  of  Uf  also  affects  the  prediction  of 
slant  fracture  during  Mode  I  loading.  Assuming  Uf  = 
le  —  5  mm  leads  to  a  prediction  of  flat  fracture  even 
though  the  mesh  size  used  is  very  fine  (see  Fig.  17a).  At 
the  same  time,  slant  fracture  is  predicted  when  using  a 
larger  value  Uf  =  0.1  mm  (Fig.  17b),  which  is  physi¬ 
cally  more  justified.  It  is  thought  that  Mode  I  response 
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Fig.  14  Influence  of  mesh 
size  on  load-displacement 
curves,  a  Mode  I, 
b  Mode  III 


Fig.  15  Influence  of  mesh 
size  on  fracture  surface  for 
Mode  I  (both  a,  b)  use 
Uf  =  0.1).  a  Mesh  size  = 
0.25  mmm,  b  Mesh  size  = 
0.0625  mm 


Fig.  16  Influence  of  Uf  on 
load-displacement  and 
fracture  surface 


Fig.  17  Influence  of  Uf  on 
fracture  surface  for  Mode  I 
loading  (both  a,  b)  use  fine 
mesh  size,  0.00625  mm), 
a  Flat  fracture, 

Uf  =  le-5  mm,  b  Slant 
fracture,  Uf  =  0.1  mm 
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Fig.  18  Crack  extension  vs.  actuator  displacement  predicted  by 
FE  simulations 

provide  a  better  environment  for  the  void  growth  and 
linkage  which  may  take  longer  to  lead  to  the  final  sep¬ 
aration.  As  a  result,  the  parameter  u  /  which  controls 
the  post-initiation  range,  should  be  larger.  By  contrast, 
Mode  III  loading  promotes  shear  localization  mecha¬ 
nism,  which  occurs  more  suddenly.  Similar  trend  was 
observed  by  Li  and  Wierzbicki  (2010a)  when  simulat¬ 
ing  plane  strain  fracture  of  TRIP  steel. 

Based  on  these  studies,  it  appears  that  the  mag¬ 
nitude  of  the  softening  parameter  should  be  adjusted 
depending  on  the  magnitude  of  stress  triaxiality.  In 
other  words,  the  rate  of  element  softening  depends  on 
the  stress  state  at  the  crack  tip.  One  can  give  here  an 
analogy  with  the  classical  concept  of  Rice  and  Tracey 
(1969),  where  equivalent  strain  to  fracture  was  found 
to  be  the  diminishing  function  of  stress  triaxiality  (high 
negative  hydrostatic  pressure) 

11  Crack  propagation 

11.1  Crack  length 

One  of  the  parameters  that  could  be  measured  from 
the  test  is  the  amount  of  surface  crack  extension,  see 
for  example,  Simonsen  and  Tornqvist  (2004).  Figure  1 8 
presents  the  predicted  crack  length  as  a  function  of  actu¬ 
ator  displacement  for  all  four  loading  angles.  The  pre¬ 
dicted  crack  extension  is  almost  linear.  It  is  interesting 
to  note  that  the  actuator  displacement-crack  extension 
results  for  60°  and  90°  almost  coincide.  Also,  the  tra¬ 
jectories  for  0°  and  30°  are  very  similar.  At  the  same 
time,  the  load-displacement  responses  for  all  four  cases 
are  quite  distinct. 


Fig.  19  Evolution  of  triaxiality  and  equivalent  strain  to  fracture 
for  Mode  I  and  Mode  III 


11.2  Evolution  of  stress  triaxiality 

According  to  the  present  fracture  model,  ductile  frac¬ 
ture  is  a  strong  function  of  stress  triaxiality  Y)  and  Lode 
angle  0 ,  which  control  the  amount  of  stress  constraint  in 
the  crack  tip  region.  For  the  present  plane-stress  prob¬ 
lem,  the  Lode  angle  is  uniquely  related  to  the  stress 
triaxiality  (see  Bai  and  Wierzbicki  2010  and  Li  and 
Wierzbicki  2010a  for  details  regarding  this  relation¬ 
ship),  so  it  is  sufficient  to  study  the  evolution  of  one 
parameter,  which  is  taken  to  be  ?].  The  stress  triaxi¬ 
ality,  ?],  is  a  function  of  spatial  position  and  loading 
history.  The  history  of  stress  triaxiality  for  an  element 
located  at  the  intersection  of  the  symmetry  planes  and 
the  pre-crack  tip  is  shown  in  Fig.  19.  The  trajectory 
lines  approach  the  fracture  locus  defined  by  Eq.  3,  and 
then  cross  it  before  fracture  initiates.  In  Mode  I  load¬ 
ing,  the  stress  triaxiality  is  approximately  constant  but 
never  reaches  the  value  0.577  corresponding  to  plane 
strain  (the  pointed  black  dash  line  in  Fig.  19).  At  the 
same  time,  for  Mode  III  loading  the  parameter?]  starts 
from  zero  (as  expected  for  pure  shear  loading)  and  then 
increases  to  a  final  value  of  0.414  at  the  onset  of  frac¬ 
ture.  This  indicates  that  there  is  a  gradual  transition 
from  pure  shear  at  the  beginning  of  the  loading  process 
to  combine  shear  and  tension  later  on.  In  other  words, 
it  is  impossible  to  generate  pure  Mode  III  loading  in 
the  present  type  of  tearing  test,  because  there  is  a  grad¬ 
ual  transition  from  Mode  III  to  combined  Mode  I/III 
loading. 
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Fig.  20  Evolution  of 
triaxiality  for  Mode  I  and 
Mode  III  at  four  different 
locations  along  the  crack 
path,  a  Mode  I,  b  Mode  III 


Actuator  displacement  (mm)  Actuator  displacement  (mm) 


Actuator  displacement  (mm) 


Fig.  21  Snapshots  of  cross  sections  at  four  locations  showing 
necking,  tunneling  and  shear  band  formation  at  the  tip  of  the 
propagating  crack  for  Mode  I  loading 


mate  at  0, 1/4, 1/2  and  3/4  width  between  the  pre-crack 
tip  and  the  end  of  the  specimen.  For  Mode  I,  approx¬ 
imately  all  elements  undergo  similar  histories  of  the 
parameter  rj.  A  sharp  peak  corresponds  to  the  post¬ 
initiation  behavior  after  which  the  equivalent  stress  a 
diminishes  from  the  maximum  value  to  zero  while  the 
mean  stress  is  approximately  constant.  Then,  accord¬ 
ing  to  the  definition  of  stress  triaxiality,  the  parameter 
t)  oo  and  then  the  element  is  deleted  immediately.  A 
different  situation  is  observed  for  Mode  III.  Elements 
close  to  the  crack  tip  starts  to  deform  under  pure  shear. 
Subsequently,  twisting  of  the  specimen  takes  place  and 
the  element  located  further  away  from  the  initial  crack 
tip  undergoes  compression  before  returning  back  to 
tension,  see  Fig.  20b. 


The  evolution  of  stress  triaxiality  for  elements 
located  along  the  crack  path  is  shown  in  Fig.  20.  Four 
elements  are  chosen  for  the  analysis,  located  approxi- 


11.3  Detailed  analysis  of  specimen  deformation  at  the 
crack  front 

The  results  of  the  numerical  simulation  provide  addi¬ 
tional  insight  into  the  details  of  the  predicted  crack 


Fig.  22  Details  of  the 
FE-captured  behavior  at  the 
crack  front  for  Mode  I 


Tunneling 


Two  types  of  necking 
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Actuator  displacement  (mm) 


Fig.  23  Twists  of  the  specimen  in  Mode  III  simulation 


formation  and  propagation  process.  Snapshots  of  the 
deformed  cross  sections  at  four  locations  (the  same 
locations  as  those  shown  in  Fig.  20)  along  the  width  of 
the  specimen  for  Mode  I  loading  are  shown  in  Fig.  21. 
From  the  deformed  finite  element  mesh,  one  can  clearly 
see  a  gradual  transition  from  flat  to  slant  fracture,  for¬ 
mation  of  a  shear  band,  local  neck  formation  ahead  of 
the  crack,  and  crack  tunneling.  Those  features  can  be 
detected  when  examining  a  partially  fractured  speci¬ 
men. 

A  summary  of  various  features  of  crack  propaga¬ 
tion  captured  by  the  FE  simulations  is  given  in  Fig.  22, 
including  necking,  tunneling  and  flat  to  slant  transi¬ 
tion.  Similar  features  were  also  predicted  by  Xue  and 
Wierzbicki  (2009)  in  the  numerical  simulation  of  CT 
specimen  using  fracture  coupled  with  damage  plastic¬ 
ity  approach. 

Different  from  Mode  I,  little  thickness  reduction  is 
observed  for  Mode  III  loading,  as  shown  in  Fig.  23. 
The  elements  around  the  crack  tip  are  deleted  in  the 
horizontal  direction  along  the  maximum  shear  starting 
from  the  free  surface.  At  the  same  time,  there  is  a  grad¬ 
ual  evolution  of  the  twist  of  the  specimen  which  causes 
a  transition  from  a  nearly  pure  Mode  III  into  a  com¬ 
bined  Mode  III/I.  It  can  be  concluded  that  the  loading 
fixtures  and  thin  sheet  material  will  not  readily  generate 
pure  Mode  III  conditions  but  rather  a  combination  of 
two  modes.  This  conclusion  is  consistent  with  a  much 
earlier  observation  of  Zhou  and  Wierzbicki  (1996)  that 
shear  punching  can  be  regarded  as  tensile  deformation 
in  the  rotated  coordinate  system. 

The  final  deformation  predicted  by  FE  simulation  of 
the  specimens  under  Mode  III  loads  is  compared  with 


the  photographs  of  the  deformed  specimen  in  Fig.  24, 
from  different  views.  Because  of  the  finite  stiffness  of 
the  gripping  system  ( ),  the  initially  parallel  edges 
of  the  specimen  are  subjected  to  twisting.  The  amount 
of  calculated  final  twist  compares  very  well  with  the 
actual  specimen. 

A  comparison  of  the  evolution  of  damage  indicator 
through  the  thickness  for  both  Mode  I  and  Mode  III 
before  fracture  at  the  cross-section  in  the  middle  of  the 
specimen  is  shown  in  Fig.  25.  For  Mode  I,  damage  ini¬ 
tiates  at  the  center  and  propagates  outward  through  the 
thickness,  whereas  for  Mode  III,  damage  initiates  from 
the  surface  and  then  propagates  towards  the  center. 

12  Discussion  and  conclusions 

The  quasi- static  crack  propagation  in  a  modified  CT 
specimen  under  combined  Mode  I  and  Mode  III  loads  is 
studied  both  experimentally  and  computationally.  The 
MMC  fracture  model  is  used  to  define  a  fracture  initi¬ 
ation  threshold  for  each  element  and  element  deletion 
is  used  in  conjunction  with  a  post-initiation  softening 
model  to  predict  crack  growth  along  the  experimen¬ 
tally  observed  crack  path.  It  should  be  noted  that  the 
MMC  model  is  not  an  ad  hoc  postulated  model  but 
a  result  of  carefully  planned  and  executed  tests  along 
with  a  suitable  theoretical  development.  Furthermore, 
this  model  is  now  widely  used  in  the  automotive  and 
steel  industry  for  predicting  initiation  and  propagation 
of  cracks  and  sheets.  The  results  of  the  present  study 
indicate  that  by  using  the  MMC  fracture  model  with 
appropriate  modifications  of  the  gripping  conditions, 
the  load  displacement  response  for  all  combined  load¬ 
ing  cases  (0°,  30°,  60°  and  90°)  can  be  simulated  with 
good  accuracy.  Also,  all  the  local  features  of  the  speci¬ 
men  deformation  around  the  crack  tip  such  as  thinning, 
tunneling  and  transition  from  flat  to  slant  fracture  were 
predicted  with  a  great  realism. 

It  was  shown  that  slight  variations  in  boundary  con¬ 
ditions  will  have  a  major  effect  on  the  transfer  of  load 
from  the  grips  to  the  crack  tip  region,  and  hence  on  the 
load-displacement  response  of  the  specimen.  Specifi¬ 
cally,  it  is  found  that  the  load-displacement  response  of 
Mode  I  is  very  sensitive  to  the  planar  rotational  degrees 
of  freedom  whereas  the  Mode  III  loading  configuration 
exhibits  less  sensitivity  to  boundary  condition  varia¬ 
tions,  though  the  effects  need  to  be  taken  into  consid¬ 
eration  when  modeling  the  entire  specimen  and  grip 
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Fig.  24  Comparison  of 
specimen  deformation  of  FE 
simulations  and 
Experiments  of  Mode  III 
a  simulation,  b  experiment 


combination.  Given  the  complexity  of  any  joint  config¬ 
uration,  the  simulations  suggest  that  a  hybrid  approach 
involving  (a)  a  combination  of  experimental  measure¬ 
ments  (3D  displacements  on  both  sides)  on  the  spec¬ 
imen  in  a  region  removed  from  the  crack  tip  and  (b) 
finite  element  analysis  using  smoothed  displacement 
data  input  on  a  pre-defined  boundary  contour  enclosing 
the  crack  tip  and  an  assumed,  through-thickness  varia¬ 
tion  (e.g.,  linear  variation)  offers  a  way  to  mitigate  the 
need  to  model  the  complex  gripping  conditions.  At  the 
same  time,  the  boundary  conditions  provide  sufficient 
input  data  to  determine  crack  tip  conditions  with  high 
fidelity.  A  similar  approach  was  recently  demonstrated 
by  Yan  et  al.  (2009a, b). 

An  extensive  parametric  study  demonstrated  that  (a) 
solid  elements  must  be  used  for  fracture  simulations  in 
order  to  capture  such  effects  as  slant  fracture,  neck¬ 
ing,  tunneling,  local  twisting  of  the  specimen  under 
out-of-plane  loading  etc.,  (b)  mesh  size,  element  aspect 
ratio  and  post-initiation  parameter  Uf  must  be  chosen 
carefully  to  predict  formation  of  slant  fracture.  In  par¬ 
ticular,  the  magnitude  of  the  post-initiation  parameter 
u f  must  be  adjusted  to  account  for  the  anticipated  active 
ductile  failure  mechanism  (void  growth  and  linkage  of 
voids  or  shear  failure),  suggesting  that  uf  is  not  a  unique 
parameter  but  rather  is  a  function  of  local  conditions 
such  as  stress  triaxiality. 

Finally,  our  simulations  demonstrate  that  it  is  not 
possible  to  generate  a  pure  Mode  III  response  using 


Mode  III 


Fig.  25  Contour  of  the  damage  evolution  through  thickness  for 
Mode  I  and  Mode  III 


the  present  type  of  specimen  and  fixture.  During  the 
tearing  process  of  thin  sheets,  the  material  is  subject 
to  torsional  loading  result  in  crack  conditions  that  are 
a  combination  of  modes,  including  both  Mode  I  and 
Mode  III. 
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A  Appendix 

A.l  Estimation  of  the  stiffness  K\  with  lower  grip 
rotation  in  clevis 

All  elements  connecting  the  specimen  to  the  hydraulic 
grips  are  relatively  solid  and  thus  large  bending  rigidity. 
The  rotational  spring  stiffness  K\  can  be  calculated  by 
considering  a  simple  model  of  a  rigid  fixture  pressing 
onto  one  side  of  an  elastic  short  cantilever  beam,  see 
Fig.  26. 

During  bending,  the  rigid  fixture  is  exerting  a  force  P 
at  the  tip  of  the  clevis.  A  part  of  the  loading  is  also  trans¬ 
mitted  by  the  pin.  Consideration  of  the  elastic  defor¬ 
mation  of  the  pin  will  lead  to  a  statically  indetermined 
problem  which  could  be  only  solved  numerically.  In 
the  simple  derivation,  it  is  assumed  that  the  only  con¬ 
tribution  of  the  bending  moment  M  at  the  root  of  the 
short  cantilever  comes  from  the  contact  force  P : 


bh3 

~12' 


(A3) 


Thus,  the  bending  stiffness  of  the  clevis  can  be  writ¬ 
ten  in  the  form; 


M 

K1  =  J  (A4> 

Combining  Eqs.  (A1)-(A4),  for  the  case  without  a 
pin, 


K 1 


Ebh 3 
61 


(A5) 


For  the  clevis  used  in  the  experiments,  the  dimen¬ 
sions  are  /  =  b  =  76.2  mm  and  h  =  2.5  mm.  Since 
Young’s  modulus  for  the  clevis  material  17-4  stain¬ 
less  steel  is  E  =  210  GPa,  the  bending  component 
of  the  stiffness  matrix  of  the  clevis  becomes  K\  ~ 
5.5  x  102 Nm. 

It  is  clear  that  the  pin  contributes  substantially  to  the 
stiffness  of  the  entire  joint.  An  order  of  magnitude  anal¬ 
ysis  indicates  that  the  above  lower  bound  of  the  stiff¬ 
ness  should  be  increased  by  the  factor  of  4.  Therefore, 
in  all  subsequent  calculations,  the  bending  stiffness  is 
assumed  to  be  K\  K  2  X  103  Nm. 


M  =  PI. 


(Al) 


From  the  simple  beam  theory,  the  rotation  angle  at 
the  tip  of  the  right  arm  of  the  fork  is 
PI 2 

e  =  — ,  (A2) 

2EI 


where  I  is  the  moment  of  inertia  of  the  rectangular 
cross  section  of  the  right  arm  of  the  fork  with  width  b 
and  thickness  h  (see  Fig.  26). 


1  1 
1  h  1 

H-N 
1  1 
1  1 
i _ 1 


Fig.  26  Simple  conceptual  model  of  the  stiffness  of  the  fork 
joint 


A.2  Estimation  of  the  stiffness  K2 


Since  the  pin  makes  the  two  arms  of  the  fork  twist 
together,  then  the  polar  moment  of  inertia  for  the  rect¬ 
angular  section,  with  the  length  ratio  b/h  larger  than  5 
for  each  arm,  is  approximately 


(A6) 


The  twist  ratio  (twist  angle  for  uniform  length) 


d(j)  _  Mt 
dx  GJ 


(A7) 


where  Mt  is  the  torque,  and  G  is  the  shear  modulus 
G  =  2(\+v)  •  Therefore,  the  torsional  component  of  the 
stiffness  of  the  fork  (including  both  arms)  is: 


K2 


2GJ  _  Eh 3 
l  3(1  +  v) 


(A8) 


For  the  present  geometry,  the  torsional  stiffness 
becomes  K2  ~  800  Nm.  It  should  be  noted  that  the 
torsional  stiffness  does  not  influence  Mode  I  loading 
response,  while  for  combined  loading  cases,  the  defor¬ 
mation  of  the  specimen  depends  on  both  K\  and  K2. 
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